Uoj 62 怎样跑得更快

Posted on | Views:

莫比乌斯反演.

为了方便,假定我们要解的方程组的形式为
$$
b_i\equiv \sum_{j=1}^n \gcd(i,j)^c {\rm lcm}(i,j)^d z_j \pmod {998244353}
$$
变形一下,得到
$$
b_i\cdot i^{-d} \equiv \sum_{j=1}^n \gcd(i,j)^{c-d} (j^d\cdot z_j) \pmod {998244353}
$$
那么方程组可以写成
$$
y_i \equiv \sum_{j=1}^n \gcd(i,j)^w x_j
$$
这样的形式,我们只要解出了这里了 $x_j$ ,由 $j^d\cdot z_j\equiv x_j$ 就可以求得 $z_j$ 了.

先尝试推一波式子把 $\gcd$ 去掉.
$$
\begin{aligned}
y_i &= \sum_{j=1}^n \gcd(i,j)^w x_j \\
&=\sum_{d|i} d^w\sum_{j=1}^{\lfloor\frac n d\rfloor}x_{jd}\sum_{s|\frac{i}{d},s|j} \mu(s)
\end{aligned}
$$
我们直接枚举 $T=ds$ ,得到
$$
y_i=\sum_{T|i}(\sum_{d|T}d^w\mu(\frac Td)) \sum_{j=1}^{\lfloor\frac n T \rfloor} x_{jT}
$$
$y_i$ 是已知的,我们用一次莫比乌斯反演即可求出每个 $T$ 的 $(\sum_{d|T}d^w\mu(\frac Td)) \sum_{j=1}^{\lfloor\frac n T \rfloor} x_{jT}$ .

前面那个式子 $\sum_{d|T}d^w\mu(\frac Td)$ 是一个关于 $T$ 的积性函数,可以线性筛预处理,也可以枚举 $d$ ,将贡献加到每个 $T$ 上.

那么两者相除就可以得出每个 $T$ 的 $\sum_{j=1}^{\lfloor\frac n T \rfloor} x_{jT}$ ,我们从大到小枚举 $T$ ,就可以得出每个 $x_i$ ,进而求出 $z_i$ .

时间复杂度 $O(n\log n)$ .

形如 $y_i=\sum_{j} f(\gcd(i,j))\cdot g(i)\cdot h(j)\cdot x_i$ 的方程组都可以类似求解. 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int read()
{
	int out = 0, sgn = 1;
	char jp = getchar();
	while (jp != '-' && (jp < '0' || jp > '9'))
		jp = getchar();
	if (jp == '-')
		sgn = -1, jp = getchar();
	while (jp >= '0' && jp <= '9')
		out = out * 10 + jp - '0', jp = getchar();
	return out * sgn;
}
const int P = 998244353;
int add(int a, int b)
{
	return a + b >= P ? a + b - P : a + b;
}
void inc(int &a, int b)
{
	a = add(a, b);
}
int mul(int a, int b)
{
	return 1LL * a * b % P;
}
int fpow(int a, int b)
{
	int res = 1;
	while (b)
	{
		if (b & 1)
			res = mul(res, a);
		a = mul(a, a);
		b >>= 1;
	}
	return res;
}
const int N = 1e5 + 10;
int c, d, w, f[N], g[N], h[N], pwd[N];
int mu[N], prime[N], cnt = 0, ism[N];
void init(int n)
{
	ism[1] = 1, mu[1] = 1;
	for (int i = 2; i <= n; ++i)
	{
		if (!ism[i])
		{
			mu[i] = P - 1;
			prime[++cnt] = i;
		}
		for (int j = 1; j <= cnt && i * prime[j] <= n; ++j)
		{
			ism[i * prime[j]] = 1;
			if (i % prime[j] == 0)
				break;
			mu[i * prime[j]] = add(P, -mu[i]);
		}
	}
	for (int d = 1; d <= n; ++d)
	{
		int pw = fpow(d, w);
		for (int T = d; T <= n; T += d)
			inc(g[T], mul(pw, mu[T / d]));
	}
	for (int i = 1; i <= n; ++i)
	{
		g[i] = fpow(g[i], P - 2);
		pwd[i] = fpow(i, P - 1 - d);
	}	
}
int n, m, x[N], y[N], z[N];
int main()
{
	n = read(), c = read() % (P - 1), d = read() % (P - 1), m = read();
	w = c - d;
	if (w < 0)
		w += P - 1;
	init(n);
	while (m--)
	{
		for (int i = 1; i <= n; ++i)
		{
			int b = read();
			y[i] = mul(b, pwd[i]);
			f[i] = h[i] = x[i] = 0;
		}
		for (int i = 1; i <= n; ++i)
			for (int j = i; j <= n; j += i)
				inc(f[j], mul(y[i], mu[j / i]));
		bool flag = true;
		for (int i = 1; i <= n && flag; ++i)
		{
			if (!g[i])
			{
				if (!f[i])
					h[i] = 0;
				else
					flag = false;
			}
			else
				h[i] = mul(f[i], g[i]);
		}	
		if (!flag)
		{
			puts("-1");
			continue;
		}
		for (int i = n; i >= 1; --i)
		{
			int sum = 0;
			for (int j = i + i; j <= n; j += i)
				inc(sum, x[j]);
			x[i] = add(h[i], P - sum);
		}
		for (int i = 1; i <= n; ++i)
			z[i] = mul(x[i], pwd[i]);
		for (int i = 1; i <= n; ++i)
			printf("%d ", z[i]);
		puts("");
	}
	return 0;
}